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Abstract 

We study the distribution of interfacial separations P{u) at the contact re- 
gion between two elastic solids with randomly rough surfaces. An analytical 
expression is derived for P{u) using Persson's theory of contact mechanics, and 
is compared to numerical solutions obtained using (a) a half-space method based 
on the Boussinesq equation, (b) a Green's function molecular dynamics tech- 
nique and (c) smart-block classical molecular dynamics. Overall, we find good 
agreement between all the different approaches. 

Keywords: Contact mechanics, randomly rough surfaces, elastic solids, 
pressure distribution, interfacial separation 



1. Introduction 

Modeling the contact mechanics between elastic solids with surfaces that 
are rough on multiple length scales is a challenging task. To perform such a 
task several theoretical approaches have been developed over the years. At the 
core of all the approaches lie approximations that relate to describing the shape 
of the contacting surfaces. In a seminal paper. Greenwood and Williamson 
( Greenwood fc WiUiamsoJ. Il966l) (GW) proposed that the contact problem 



between two elastic rough surfaces could be reduced to the problem of one 
infinitely-hard rough surface acting on a flat elastic counterface. Within their 
model, the rough topography was described by a large collection of hemispherical 
asperities of uniform radius (which individually satisfied the Hertzian approx- 
imation) with a height distribution that followed a Gaus sian law. This initial 



appro ach was later extended by Greenwood and Tripp ([Greenwood fc Tripd . 
llQTOf ) by considering the presence of roughness on the two contacting surfaces. 
Further contributions to the original GW methodology have been proposed by 
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these models rely on the definition of "asperity". The asperity concept itself 
has proved quite controversial and depends on the resoluti on of the instrument 
used to measure the surface profile (|Poon fc Bhushanl 119951) . Another drawback 
of GW-type approaches is that using only a few parameters to describe the 
surfaces generates a one-to-many mapping possibilities, i.e., the same set of pa- 
rameters can be deduced for surfaces obtained by completely different machining 
processes. 

In spite of the great increase in computing power in the past decade, ana- 
lytical theories are still very much needed to understand the contact mechanics 
of solids with surfaces that display roughness on more than three decades in 
length-scales. Theoretical models to tackle such problems rely on approxima- 
tions and idealizations in order to analytically solve the equations of elasticity. 
For example, such equations can be exactly solved for the contact problem of 
a parabolic tip acting on a flat surface under the assumptions of linearly elas- 
tic, frictionless materials (Hertz model). Similarly, an exact solution can be 
obtained for the conta ct between a sinuso idal elastic surface and a rigid plane 
(Westergaard model) ( WestergaardL Il939l ). The GW model and its extensions 
are further examples of how to deal with surface roughness in contact mechan- 
ics, resulting in simple analytic formulas. These easy-to-handle formulas have 
proved to be of great importance, and are frequently employed in the design 
process of new technical applications. The Hertz and the Westergaard mod- 
els provide accurate representation of the mechanics of single asperities. In 
addition, the GW model describes approximately the (low squeezing-pressure) 
contact between surfaces exhibiting roughness on a single or a narrow distribu- 
tion of length scales. However, most real surfaces have roughness over many 
decades of length scales. Here, the long range elastic coupling between the as- 
perity contact regions, which is negle cted in the GW model, is now known to 



strongly influence contact mechanics ( Persson , 20081 : Gampaha et al. . 2008 ). If 



an asperity is pushed downwards at a certain location, the elastic deformation 
field extends a long distance away from the asperity i nfluen cing the contact 



involving other asperities further away ([Persson et al.l . 120021 ). We note that 



the lateral coupling between contact regions is important for arbitrary small 
squeezing pressure or load. The reason for this is that surfaces with roughness 
on many length scales can be considered as consisting of large asperities popu- 
lated by smaller asperities, with the smaller asperities being populated by even 
smaller asperities and so on. Thus, when two solids are squeezed together by a 
very small external force, the distance between the macro-asperity contact re- 
gions will be very large, and one may be tempted to neglect the elastic coupling 
between the macro-asperity contact regions. However, the separation between 
the micro-asperity contact regions within a macro-asperity region will in gen- 
eral be very small, and one cannot neglect the elastic coupling between such 
micro-asperity contact regions. This latter effect is neglected in the GW theory, 
signiflcantly limiting its prediction capabilities when applied to most real sur- 
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faces. Additionally, in the GW model the asperity contact regions are assumed 
to be circular (or elliptical) while the actual contact regio ns (at high enough ex- 



perimental resolution) show fractal-like bou ndary lines (jBorri-Brunetto et al 



200lt iPei et all 120051: iPersson fc Yan^ |2008| ). Therefore, because of their com- 



plex geometries, one should try to avoid explicitly invoking the nature of the 
contact regions when searching for an analytical methodology to solve the con- 
tact problem of two elastic rough surfaces. 

Recently, an analytical contact mechanics model that does not use the asper- 
ity concept and becomes exact in the limit of complete contact has been devel- 
oped by Persson (jPerssonl l200ll l2006l 120071: lYang fc Perssoni 120081: iPersson fc Yanei 
20081) . The theory accounts for surface roughness on all relevant length scales 



and includes (in an approximate way) the long range elastic coupling between 
asperity contact regions. In this theory the information about the surface enters 
via the surface roughness power spectrum C(q), which depends on all the sur- 
face roughness wavevectors q components. The theory can be used to calculate 
the interfacial stress distribution P{aX), from which one can obtain the area 
of real contact as a function of the squeezing pressure p and the magnification 
Furthermore, the theory predicts the average interfacial separation u for any 
applied external load. 

Besides analytical approaches, numerical algorithms (deterministic) have 
also been developed to understand the contact mechanics of elastic solids with 
rough boundaries. As the speed and memory capacity of computers increase, 
numerical methods have become a viable alternative to analytical methods when 
modeling surfaces of three-dimensional (3D) solids having surface roughness ex- 
tending over at most three decades in length-scales. Nevertheless, simplifying 
assumptions about the material and the topography are still needed to ensure 
reasonable computational time windows. Much is still to be done in order to 
reach the capacity to numerically simulate real surfaces that may have rough- 
ness from the nanometer scale up to the macroscopic size of the system which 
could be cm. 

Numerous numerical works have been reported in the literature aiming to 
solve the contact mechanics of two linear elastic solids with rough surfaces. 
The majority of these are half-space models in which the elastic deformation 
is related to the stress field at the surfaces of the solids through integral equa- 
tions where the domain of integration is the boundary of the half-space. This 
type of approach is commonly referred to as boundary element method (B EM). 
Twenty years ago Lubrecht and loannides ( Lubrecht fc loannides . 199ll ) sug- 
gested applying multilevel techniques to facilitate the numerical solutio n of the 
BEM. With the same objective in mind Ren and Lee ( Ren fc Led . 1994 ) imple- 
mented a moving grid method to reduce storage of the influence matrix when 
the conventional matrix inversi on approach is used to solve this type of prob- 
lem. Bjorklund and Andersson ([Andersson fc S. Biorklundl 119941 ) extended the 
conventional matrix inversion approach by incorporating friction induced defor- 
mations. Alternative techniques that aim to solve the elastic contact of rough 
surfaces are t he Fast Four i er Tr ansform (FFT)-based method introduced by 
Ju and Farris ( Ju fc Farrisl 19961 ) and a follow-up extension, based on a varia- 
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tional principle ( Kalker . 1977 ). proposed by Stanley and Kato ( Stanley &: Katol 
)97|) . The contact between solids with realistic surface topographies under rel- 
atively small loads usua lly leads to plastic deformations. Tian and Bhushan, 
( Tian fc Bhushanl . ll996l) Isased their theoretical model is based on a variational 
principle for linear elastic perfectly plastic materials. In this way, not only the 
in-contact topography and the corresponding pressure distribution but also the 
unloaded plastically deformed topography can be obtained for the case when 
the loading is hi gh enough to cause yield. This model was further developed 
in the paper by dSahlin et al.l . l2010l) and this is also the BEM employed in the 
present work, but here we restrict the analysis to linear elastic materials. 

In earlier works, the prediction of Persson's contact mechanics theory for 
the interfacial stress distribution P{<t) and the contact area have been com- 
pared to numerica l results obtained usin g the finite elernent ni ethod (FEM) 
( Hvun et al. . 2004 ). molecular dynamics (lYan g fc Persson , l200 8l) and Green's 
function molecular dynamics (GFMD) ( Campa iia fc Mtiseil . 12007). In this pa- 
per, we will show how the theory can be extended to also predict the distribution 
of interfacial separations P (u). This quantity is of crucial importance for pro b- 
lems like leak-rate of s eals (IPersson fc Yanel. 2008 : Lorenz fc Perssonl 2010allbh 
or mix ed lubrication ( Persson . I2OIOI : Scaraggi et al. . 20 lit iLorenz fc Persson , 
2010d) . The analytical results will be compared to numerical solutions obtained 



using (a) a half-space method based on the Boussinesq equation (BEM), (b) a 
Green's function molecular dynamics technique and (c) smart-block molecular 
dynamics. 



2. Contact mechanics theory of Persson 

Consider the frictionless contact between two elastic solids with Young's 
elastic moduli Eq and Ei and Poisson ratios Vq and i^i . Assume that the surfaces 
of the two solids have height profiles /io(x) and /ii(x), respectively. The elastic 
contact mechanics for the solids can be mapped into that of a rigid substrate 
with height profile /i(x) — /io(x) + ft.i(x) and a second elastic solid with a flat 



surface and Young's modulus E and Poisson ratio v chosen so that (| Johnson , 

mi) 

i^ = i^ + i^. (1) 

E Eq El 

The main physical variables that characterize the contact between the solids 
are the stress probability distribution P{(j) and the distribution of interfacial 
separations P{u). These functions are defined as follows: 

P(a)^(5[a-a(x)]), P{u) ^ {5[u - u{^)]) 

where 5{..) is the Dirac delta function, and ct(x) and u(x) are the stress and the 
interfacial separation at point x — (x, y), respectively. The (..) brackets denote 
ensemble averaging. Note that both P{a) and P{u) have a delta function at 
the origin with its weight determined by the area of real contact i.e. given by 
(1 — A/Aq)5{(t) and {A/ Aa)5{u). Here Aq is the nominal contact area and A 
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the area of real contact projected on the xy-plane. NormaUzation conditions 
require that 



J dcrP{a) = 1, J duP{u) = 1 



while 

daP{a)^—, / duP{u)^^^. 

Thus from the interfacial distribution of stresses or separations one can immedi- 
ately determine the area of real contact A. The average interfacial stress (which 
must be equal to the applied pressure) a, and the average interfacial separation 
u, can be obtained as 



a — J da aP{a), u — j du uP{u) 



The stress and interfacial separation distribution functions, P{a) and P{u), 
are determined by the elastic energy Ud stored i n the asperity contact regi ons 
(see below). The elastic energy U^i is written as (|Persson[ [20021 l2006l [2008h 

4(1 -z/2) 



Uei = 77^^ / <lC{q)W{q) (2) 



where the surface roughness power spectrum is defined by 

Ciq) ^j^Jd'x (Mx)MO))e-'>-. (3) 

The height profile /i(x) of any rough surface can be measured routinely nowadays 
on all relevant length scales using optical and stylus experiments. 

For complete contact W(q) = 1 render i ng an exact result for the expres- 



sion of the energy above. In Ref. (jPerssonl 120021 ) it was argued that W{q) 



P{q) = A{()/Ao is the relative contact area when the interface is studied at 
the magnification C = q/qo (where qo is the small-wavevector cut-off, usually 
chosen as tt/L, where L = ^^Aq is the linear size of the surface). The qualitative 
explanation for such an argument is that the solids will mainly deform in the 
regions where they make contact, and most of the elastic energy will arise from 
the contact regions. Nevertheless, using W{q) = P{q) assumes that the energy 
(per unit area) in the asperity contact regions is just the average elastic energy 
(per unit area) as if complete contact would occur. This does not take into ac- 
count that the regions where no contact occurs are those regions where most of 
elastic energy (per unit area) would be stored if complete contact would occur. 
Hence, we expect smaller stored elastic energy (per unit are a) in the asperity 



contact regions than obtained using W(q) = P{q)- In Ref. (lYang fc: Persson . 



2008t IPerssonl l2008t ICampafia et all . 120111) we found that using 



W{q) = P{q) [7 + (1 - l)P^{q)] = P{q)S{p, g), (4) 

with 7 « 0.45 gives good agreement between theory and numerical calculations. 
Note that for complete contact P{q) = 1 and hence W{q) = 1 which reduces to 
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the exact result for the elastic energy in such a limit. On the contrary, in the 
limit of small contact, P{q) « 1 which yields W{q) ~ lP{q)- For 7 « 0.45 
this results in an elastic energy which is a factor of 0.45 smaller than the elastic 
energy (per unit area) stored in the contact region in the c ase of complete 
conta ct. Recently, in an independent study, Akarapu et al. ([Akarapu et all . 
2OIOI) found a value of 7 = 0.48 after analyzing a variety of rough surfaces in 



contact with roughness down to the atomic scale, variable Poisson ratio and 
Hurst exponents oi H = 0.5 and 0.8. 



T he contact mechanics for malism developed by Persson (|Perssonll2001ll2006 



2007t lYang fc Perssonl 120081 ) is based on studying the interface between two 
contacting solids at different magnifications When the system is studied 
at the magnification ( it appears as if the contact area (projected on the xy- 
plane) equals A{(), but when the magnification increases, it is observed that 
the contact is incomplete, and the surfaces in the apparent contact area A{() 
are in fact separated by the average distance u{(), see Fig. [5] The (apparent) 
relative contact area A (C) /Aq at the magnification C is given by ( Perssonl I2OOII : 
Yang fc Pers"so3 . 12OO8I) 



where po = Fj^ /Aq is the nominal squeezing pressure and 

2 rCqo 



E 



1 



dqq'C{q)S{p,q). 



90 



In most applications A/Aq << 1 and in this case one may use S 
distribution of interfacial stress is given by (for a > 0): 



Pia) 



1 



2(7rG)i/2 



exp 



AG 



cxp 



4G 



(5) 

(6) 
7. The 

(7) 



Let us define ui(C) to be the (average) height separating the surfaces which 
appear to come into contact when the magnification decreases from ^ to C — A^, 
where AC is a small (infinitesimal) change in the magnification. wi(C) is a 
monotonically decreasing function of and can be c alculated from the aver age 
interfacial separation u{C,) and A{C) using (see Ref. (|Yang fc Perssonll2008l) ) 



u,{C)=u{0+u'{C)A{C)/A'{C), 



where 



uiO dq q^C{q)w{q) / dp' -S{p\q)e-^ 



Aril 

E* = E/{1 - 1,^), p{C) = poAo/A{C) and 

w{qX)= r dq' q"C{q') 



w(qX)p'/E'f 



(8) 



(9) 



-1/2 
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Figure 1: An elastic block (dotted area) in adhesive contact with a rigid rough 
substrate (dashed area). The substrate has roughness on many different length 
scales and the block makes partial contact with the substrate on all length scales. 
When a contact area is studied, at low magnification it appears as if complete 
contact occurs, but when the magnification is increased it is observed that in 
reality only partial contact has taken place. 



magnification ^ 



elastic solid 




Figure 2: An asperity contact region observed at the magnification (. It ap- 
pears that complete contact occurs in the asperity contact region, but when the 
magnification is increased to the highest (atomic scale) magnification ^i, it is 
observed that the solids are actually separated by the average distance u{(). 
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Figure 3: An elastic block squeezed against a rigid rough substrate. The 
separation between the average plane of the substrate and the average plane of 
the lower surface of the block is denoted by u. Elastic energy is stored in the 
block in the vicinity of the asperity contact regions. 

The definition of the distribution of interfacial separations P{u) — {6[u — 
u(x)]) involves an ensemble average over many realizations of the surface rough- 
ness profile. If the surface roughness power spectra has a roll-off wavevector qc 
which is much larger than go — t^I L, where L is the linear size of the surface, 
then performing an ensemble average is identical to averaging over the surface 
area. In this case we can write the distribution of interfacial separations as 

^(") = ^ / '^'^^ ~ "W)- (10) 

^0 Jao 

The probability distribution is normalized 



j du P{u) = 1. (11) 



In the contact mechanics theory of Persson (jYang fc Persson , 2008h the interface 



is studied at different magnification C. As the magnification increases, new short 
length scale roughness can be detected, and the area of (apparent) contact A{() 
therefore decreases with increasing magnification. The (average) separation 
between the surfaces in the surface area which (appears) to move out of contact 
as the magnification increases from ^ to C + ^C; is denoted by ui(C) and is 
predicted by the Persson theory (see above). The contact mechanics theory of 
Persson does not directly predict P(u) but rathe r the probability distribution 
of separation ui (see Ref. ( Yang fc Perssonll2008l )): 



A(^)-^^ 'dCM'(C)] S{u-u,{C)). (12) 

Since ui(C) is already an average, the distribution function Pi{u) will be more 
narrow than P{u), but the first moment of both distributions coincide and is 
equal to the average surface separation: 



/•OO /'C 

/ du uP{u) = 
Jo Jo 



du uPi{u). 
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To derive an approximate expression for P{u) we write Eq. (|10p as 

Piu) = ^ 1^^' dCi-A'iO] {S{u - «(x)))c. (13) 

Here {..)( stands for averaging over the surface area which moves out of contact 
as the magnification increases from C, to ( + d(^. Note that 

(u(x))c-ui(C). (14) 

A surface which moves out of contact as the magnification increases from C 
to C + will have short-wavelength roughness with wavevectors larger than 
q > C<?o- Thus the separation between these surface areas will not be exactly 
'Ui(C), but will fluctuate around this value. One may take this into account by 
using 

((u(x)-ui(C))2)c«/»L.(C), (15) 
where /irins(C) is the mean of the square of the surface roughness amplitude 
including only roughness components with the wavevector q > qqC- We can 
write 

hLsiC) = / d'l C{q), (16) 

J q>qoC 

where the surface roughness power spectra C{q) can be calculated from the 
measured surface topography. Using the definition 

Ztt 

one can rewrite Eq. (IT^ as 

I dC, [-A'(C)]-!- / da e''"("-"i(C))(e^"("i(C))-"W)).. 
Aq J 27r J 

To second order in the cummulant expansion 

P^l- [ dC h^'(C)]— / da e-("-"i(C))-"^((«i(C)-«(x))^),/2 



or using Eq. (fT5|) : 



1.7 ^' '^"(2,I.L.(0)"' V 2'.?„.«) 

The above expression does not satisfy the normalization condition Eq. (ITT]) . We 
will therefore use instead 



P.i_/.CM'«)1 



1 



(2^/^?„.s(C))'/' 



(17) 



9 



The added term in this expression can be considered as resulting from the 
cummulant expansion of 



1 



(fx S{u + m(x)). 



Note that such a term vanishes for u > 0. 

The theory described above predicts that for smah squeezing pressures p 
the area of real contact is proportional to the squeezing pressure, while the 
interfacial separation depends logarithmically on p. Both results are related to 
the fact that when increasing p existing contact areas grow and new contact 
areas are formed in such a way that, in the thermodynamic limit (infinitely- 
large system), the interfacial stress distribution, and also the size distribution 
of contact spots, are independent of the squeezing pres sure as long as thes e 
distributions are normalize d to t he real contact area A |Persson et alj . 120051 ). 
In Rcf. (Lorcnz & Pcrsson, l2009f) (see also (Lor enz et all I2OIOI )') experimental 
results were presented to test the dependence of u on p. In the study a rubber 
block was squeezed against an asphalt road surface, and good agreement was 
found between the theory and experiments. The fact that A ^ p for small load 
is also well tested experimentally, and is usually considered as the explanation 
for Coulomb's friction law which states that the friction force is proportional to 
the load or normal force. 





12 5 10 20 60 100 



Figure 4: Graphical representation of a rough topography with Hurst expo- 
nent H = 0.3 and its corresponding height-height correlation function C{q) = 
(|/i(q)P) in Fourier space. The surface topography was created using a Fourier 
filtering technique and a hard cutoff qc = 64 was imposed to it in Fourier space 
(in units of 27r/L, where L is the linear size of the simulation cell). The contin- 
uous line represents the ideal algebraic scaling expected for the height-height 
correlation function of such a surface. 



3. Numerical methods 

When two elastic solids with rough surfaces come into contact, the elastic 
deformations perpendicular to the contacting plane extend into the solids a 
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characteristic length A that could be as large as the contacting plane's lateral 
size L. Thus, in order to properly capture the mechanical response of the 
solids within the contact region, the elastic properties of the material have to 
be considered up to a distance L in the normal direction to the contacting 
plane. This is why standard algorithms that use a full representation of the 
system display a computational effort that scales with the system's size L^. 
Such a scaling rapidly becomes a limitation when the linear dimension of the 
solids increases. The previous arguments explain why coarse-grained numerical 
techniques are needed when studying the contact mechanics of solids with more 
than two decades in surface roughness length-scales. 

The theory developed in Sec. [5] for the distribution of interfacial separations 
P{u) will be compared to the predictions of the following three different coarse- 
grained numerical methods: (a) a half-space method based on the Boussinesq 
equation (BEM), (b) a Green's function molecular dynamics technique (GFMD) 
and (c) smart-block classical molecular dynamics (MD). For this, we have con- 
sidered the contact between an elastic block with a flat bottom surface and a 
randomly rough rigid substrate. Self afhne fractal topographies with Hurst ex- 
ponent values oi H = 0.3, 0.5 and 0.8 (corresponding to the fractal dimension 
D{ = 3 — H = 2.7, 2.5 and 2.2) are used to model the rigid substrate. Fig. 2] 
shows the surface topography /i(x) and the power spectrum C(q) (on a log-log 
scale) of one of our surfaces with H = 0.3, hard cutoff qc = 64 (in units of 2tt/L, 
where L is the linear size of the simulation cell) and rms-slope ^ 0.03. 

Randomly rough substrate profiles were generated on a two-dimensional 
square grid with 2048 x 2048 mesh points. The surface heights were obtained 
via a Fourier Filtering Algorithm. In the case of the GFMD and BEM methods 
the elastic interactions within the original elastic block are chosen such that 
both Lame coefficients satisfy A = /i = 1. This choice of the coefficients results 
in a Young modulus of £' = 5/2, a bulk modulus oi K = 5/3, and a Poisson 
ratio of = 1/4. 

In the smart-block molecular dynamics simulations we used a smaller system 
size than in the other two numerical schemes. The surfaces were obtained by 
choosing every 4'th grid point from the original surfaces with 2048 x 2048 reso- 
lution. This procedure yielded surfaces with x Ny — 512 x 512 mesh points. 
Since the original surfaces were quite smooth at short length scales (see the 
power spectrum in Fig. the surfaces used in the smart-block MD simulations 
are expected to give the same contact mechanics as the original topographies. 
This was confirmed by comparing the MD results obtained for the 512 x 512 
system sizes to those of the BEM method for the equivalent 2048 x 2048 systems. 
In our MD simulations the atoms in the bottom layer of the block are located 
on a simple square lattice with lattice constant a — 2.6 A. The lateral dimen- 
sions of the block and substrate are L^. = Ly = N^a = 1331.2 A. The Young 
modulus of the block is E = 250 GPa and its Poisson ratio v = 1/A identical 
to that of the other two schemes. Since no natural length scale exists in elastic 
continuum mechanics, one can directly compare the results of the smart-block 
MD model to those of the methods (a) and (b) by simply using a distance scal- 
ing factor of (512/2048)a = 0.65 A, and a pressure (or stress) scaling factor of 
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250/2.5 = 100 GPa. 

Next, we will provide a short technical review of the three numerical methods 
that we have employed in this work. 



3.1. Boundary Element Method 

Any contact mechanics problem can be solved by using a technique that 
minimizes the total potential energy of the system. Assuming frictionless linear 
elastic contact, the variational pro blem includin g constraints to be solved can 
be expressed as Eg . ([T5| (see, e.g., ( Kalker , 1977 ). ( Tian fc Bhushan . 19961) and 



( Sahlin et all . 120101) 1: 



mm 

p>0 



- / d^x p(x)m2(x) - / d\ p(x)u*(x) 



(18) 



where p(x) is the pressure distribution, Uz{Ti) is the associated elastic deforma- 
tion and Wz(x) is the prescribed surface displacements, equivalent to the rough- 
ness height coordinate {h) plus a constant controlling the prescribed shift. The 
first term describes the internal complementary energy due to elastic deflection, 
and the second term governs the contribution from the prescribed displacement 
M*(x). Note that u(x) — Uz{x) — u*(x). The Boussinesq relation between 
pressure and elastic deformation employed for this work may be formulated as 



Wz(x) 



ttE 



2„/ P(x') 



(19) 



For the BEM method employed here, the complete system of equations consists 
of Eq. ([T5)) and Eq. ([T^ , and the force balance relation 



d2a;p(x)=FN, p(x)>0, 



(20) 



where -Fn is the normal force or the applied load. 

Numerica l ly, the solution is achieved by employing the method outlined by 
(|Sahlin et al.l . l2010l) . where the EFT algorithm is utilized to accelerate the com- 
putation of the elastic deflection. Eor the numerical results presented in this 
work, convergence of the solution process is reached when the following two 
convergence criteria are met: 

- The force balance criterion that controls that the load generated by the 
contact pressure supports the applied load. 



1 



d^x p(x) 



< 10" 



A geometric criterion that controls that points 'in-contact' lie sufficiently 
close to the contact plane. 



max |m2(x) 



max /i(x) 



min /i(x) 
xeAo 



< 10" 
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3.2. Greens Function Molecular Dynamics 

Greens Function Molecular Dynamics (GFMD) ( Campana fc Miiser . 20061 ) is 



among the many simulation techniques available to find the equilibrium configu- 
ration of a mechanical system under the action of external loads. To achieve its 
goal GFMD solves the system's equations of motion for a set of initial/boundary 
conditions. 

The potential energy of a linear elastic (harmonic approximation) solid is 
given by 

^ = ^ E E fc.f "»%73 = ^ E • • (21) 

ij a/3 ij 

where — UiaQa (where ei = x, 62 = y and 63 = z are orthogonal unit 
base vectors) is the deformation field at locations x^, and Kij = X]q^ k'^fsaep 
the force constant matrix. In thermal equilibrium, the displacements will 
comply with the Boltzmann statistics with second moments given by 

(uiU,) = / rfui • • • dujv Ui u,- e"''^ 



Z „ 

- kBT[K-^]^^, (22) 

with Z being the partition function. The equation above shows that one can 
obtain the force constant matrix from measuring the fluctuations in 

the second moments (ujUj). Furthermore, if the bottom surface of the solid is 
exposed to an external force field Fext({u}), all the degrees of freedom that 
do not belong to the bottom surface can be integrated out (eliminated) yielding 
an equivalent problem for the bottom surface 

surface surface 

V= 2"'-^yu.- E F-fU^ (23) 

ij i 

where Kij are new renormalized force constants. Renormalization takes place 
such that the new 2D surface obtained will deform under the action of the 
external field in exactly the same way as the bottom surface of our original 
3D-solid. The matrix 



- ij 

is known as the Greens' function of the system. From Eq. (I23p one obtains the 
equilibrium condition 

J2K^rn,=F^l (24) 

j 

If the system is periodic in the (x, j/)-plane (translational symmetry) one can 
use the Fourier transform to obtain decoupling of the modes in the q — {qx, Qy)- 
space as 

u(q) = G(q)Fext(q). (25) 
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Eq. eliminates the non-local nature of the real-space solution of Eq. (|24p 
while rendering an easily parallelizable scheme that can be used to simulate the 
contact mechanics of large systems. For a given interaction kernel G'(q), it is 
the implementation of Eq. (j25p in a molecular dynamics fashion that lies at the 
core of GFMD. 

Ou r GFMD implemen t ation followed the approac h described in previous 



works (jCampana fc Miiseii . l2007HCampana et al 



20081 ) where all the roughness 



is placed on the rigid substrate and the elasticity on a flat GFMD block. The 
interactions between the block and the rough substrate are modeled via a hard- 
wall potential. If the z-coordinate of an atom within the elastic block at location 
X = {x,y) crosses through h{x), the corresponding interaction energy increases 
from zero to infinity. To obtain the value of the surface height /i(x) inside any 
square element of the grid, we employed interpolation via bi-cubic splines with 
zero partial and cross-derivatives at the corner grid points of each element. 

Defining which block's atoms are in contact after equilibrium has been 
reached is done by analyzing the pressure distribution. In a fully equilibrated 
simulation, a wide pressure gap of several orders of magnitude would exist be- 
tween the atoms that belong to the contact region and those which do not. 
This approach to defining the contacting status of a certain atom differs from 
geometrical ones where contact is defined by comparing the relative distance be- 
tween the block's atom and the corresponding surface height at the local atomic 
position, as it is done within the BEM method. 

3.3. Smart-block Molecular Dynamics 

The smart-block molecular dynamics (MD) system is composed of an elastic 
block interacting with a rigid randomly rough substrate. The substrate and 
the bottom layer of the block consist of an array of 512 x 512 atoms. Periodic 
boundary conditions are applied in the xy-plane. The atoms in the bottom layer 
of the block form a simple square lattice with lattice constant a = 2.6 A. The 
mass of a block atom is 197 amu, and its elastic parameters have already been 
mentioned at the beginning of this section. 

In order to allow for a correct description of the long- wavelength components 
in the deformation field of the block, its thickness is chosen to be 1350.7 A, which 
is slightly larger than its lateral dimension. The techni cal details of the s mart- 
block implementation has been discussed elsewhere, see ( Yang et al. . 20061 ). The 



current smart-block consists of 12 atomic layers, and merging factors of 2 (in 
all 3 directions) are used for all layers, except the I'st, 6'th and the 11 'th. The 
smart-block contains 615780 atoms, and the total number of atoms involved in 
the simulations is 877924. 

The atoms at the block-substrate interface interact via a repulsive potential 
U{r) = 4e (ro/r)^^, where r is the interatomic distance and the parameter e cor- 
responds to the binding energy between two atoms at the separation r = 2^/^ro. 
In our calculations we have used the values tq = 3.28 A and e = 18.6 meV. Zero 
tempera ture is ma i ntaine d during the simulations using a Langevin thermo- 



stat (jGriebel et al.l . 120071 ) and the equations of motion have been integrated 
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20071: iRapapord \2004 ) with a time step 



using Verlet's method (jGriebel et al 
of At = 1 fs. 

In the present study the squeezing process proceeds as follows. The upper 
surface of the smart-block is moved towards the substrate at a constant velocity 
oi V = 5 Ta/s with the block being compressed as its bottom layer approaches 
the substrate. The duration of simulations depends on the type of substrate 
and they last until a small enough separation between the bottom block layer 
and the substrate is achieved. Note that the thickness of the smart-block may 
influence the results. In particular, a too thin smart-block leads to noisy data 
with considerable deviation from the results of the other methods. Increasing 
the thickness of the smart-block beyond the lateral size of the surface ensures 
convergence in the results. One must add that the time dependencies of p and u 
are not monotonic, and some oscillations are observed, which may be attributed 
to elastic waves propagating in the block during its compression. Lowering the 
velocity of movement down to 2.5 m/s leads to a decrease in the amplitude of 
these oscillations. Moreover, the large u region in the dependence of pressure 
on u gets closer to the BEM results when a lower value of v is employed. This 
suggests that smaller values of v should be used in future studies. 



4. Results and discussion 

We first consider the dependency p{u) of the pressure p on the average 
interfacial separation u. Figure [5] shows p{u) obtained using the analytical 
theory (black lines), the BEM (blue symbols), GFMD (green symbols) and the 
smart-block MD approach (red lines). The figure includes results for surfaces 
with Hurst exponents oi H = 0.3, 0.5 and 0.8. Note that the analytical theory 
has been developed for infinite systems, which will have infinitely high asperities, 
therefore always leading to a certain degree of contact between the solids. The 
small systems sizes utilized in the numerical simulations resulted in the highest 
asperities exceeding by only a factor of three (above the average plane) the 
root-mean-square roughness of the surfaces. This explains the sharp drop in 
the pressure in the computer simulation curves at the height threshold value 
established by the tallest asperity. 

In the smart-block MD simulations, the finite range of the interaction poten- 
tial between the block atoms and the surface atoms resulted in a non-unique way 
of defining u. In the current work, to obtain u a small contribution 6u of about 
4 A has been subtracted from the difference zi — zo in the average position of 
the interfacial atoms of the block and the substrate. How to find the best 6u is 
a rather difficult question, and several alternative ways exist to accomplish such 
a goal. Here we have used that in the continuum limit the P(u) distribution 
must display a maximum (delta-like behaviour) at its origin. In smart-block 
simulations the maximum in P{u) gets shifted to a non-zero surface separation 
due to the finite range of the interaction potential. This behavior is shown in 
Fig. [6] for two different external pressure values. Thus, shifting the MD proba- 
bility distribution P{u) towards the origin by Su is necessary in order to be able 
to compare with the GFMD and BEM continuum mechanics results. 
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Persson's contact mechanics theory predicts that the average interfacial 
separation in the large u range is related to the applied pressure p via p ~ 
exp(— u/uq), where uq is of order of the root-mean-square roughness of the 

undefo rmed rough surfac e. T his result, already confirmed by experimental 

works ( Lorenz fc PerssonLlioollLorenz et al.l . l2010l ). differs drastically from the 
prediction of GW-like asperity models which instead yield a dependency of the 
type p ~ exp(— with b being a constant. The origin for those differences, an 
exponential decay predicted by Persson's theory and a Gaussian decay obtained 
by GW, is the omission of the long-range elastic deformations within asperity 
contact models which also results in very different morphologies for the contact 
regions as illustrated in Figs. [7] and HI 

Figure [7| shows the contact morphology generated for a H — 0.8 rough 
surface at three different external pressures, p = 0.0032i?, O.OOSi?, 0.012i?, 
when the long-range elastic deformations are taken into account. In Figure[8]we 
show the morphologies obtained (for identical values of the "true" contact area A 
) when a bearing area model is utilized. The lack of long-range elastic coupling in 
the latter model produces qualitatively different contact morphologies (compare 
Figures [7] and [8]) . Thus, when elastic deformations are considered, the contact 
regions become less compacted, with fractal-like boundaries, and are distributed 
over a larger fraction of the nominal contact area than those predicted by the 
asperity model in Fig. [H 
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Figure 6: Distribution of interfacial separations for the substrate with H = 0.5 
obtained in classical MD simulations. 

A first comparison between the probability distribution P{u) of interfacial 
separations u obtained using the BEM, smart-block MD and GFMD method- 
ologies at the squeezing pressures p = 0.79, 0.78 and 0.75 GPa, respectively, is 
displayed in Figure [Q] As shown, slight changes in applied pressure did not vary 
significantly the general shape of P{u). This is indicative of the individual con- 
sistency in the implementation of each numerical technique. Further comparison 
of the numerical methods to the theoretical predictions is included in Fig. [10] 
in the limit of (a) non-contact, (b) low-pressure region and (c) medium-to-high 
pressure region. The red lines correspond to the predictions of Persson's con- 
tact mechanics theory (Eq. (|17l) ) while the blue lines have been obtained from 
GFMD (in (b)) and BEM ((a) and (c)). All the numerical results correspond to 
a single realization of the rough surface which explains the rather large fiuctua- 
tions (noise) in the P{u) data. In particular, the ensemble averaged P{u) for the 
non-contact case (zero squeezing pressure, p = 0) must follow a Gaussian law 
as given by the theory curve. However, the lack of a small wavevector roll-off 
(or cutoff) in the surface roughness power spectra implied that numerous inde- 
pendent topographies would have to be considered in order to achieve averages 
that closely represent the thermodynamic limit. Due to the large computational 
effort involved in such a task no further attempt to improve our statistics was 
performed. 

As already mentioned in the methods section, the original contact mechanics 
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(a) p/E = 0.0032. (b) p/E = 0.0080. (c) p/E = 0.0120. 



Figure 7: Contact morphologies for the H = 0.8 surface at three different 
contact pressures; (a) p = 0.0032£;, (b) p = 0.0080£: and (c) p = 0.0120£: when 
long-range elastic deformations are considered. The fractional contact areas are 
A/Ao = 0.066, 0.161 and 0.238, respectively 




(a) A/Ao = 0.066. 



(b) A/Ao = 0.161. 



(c) A/Ao = 0.238. 



Figure 8: Contact morphologies predicted by a bearing area model for the same 
"true" contact area A values as in Figure[7]but where long-range elastic coupling 
has been neglected. 
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0.12 




Figure 9: Probability distribution P{u) of interfacial separations u, obtained 
using the BEM, MD and GFMD methods for p = 0.79, 0.78 and 0.75 GPa, 
respectively. 



theory of Persson does not directly predict P{u) but instead Pi{u) (distribu- 
tion of boundary-line averaged interfacial separations) which is a much sharper 
function than the interfacial separation distribution. The individual behaviour 
of both functions for the H = 0.8 surface at a squeezing pressure oi p = 0.003E' 
is plotted in Fig. [11] In the figure the analytical P{u) predicted using Eq. (|17p 
(red curve) and Pi{u) (green line) (Eq. (|12|) ) are compared to the numerically 
exact distribution generated from GFMD simulations (blue curve) . Because the 
bin size of the P{u) and Pi{u) computations was smaller than that of the nu- 
merical GFMD study the delta function at the origin u — barely shows in the 
first two cases. Nevertheless, based on the results presented in Figs. HI [TOland 
[TT] we feel confident to conclude that the extension presented in this work to 
compute P{u) within the framework of Persson's original contact theory yields 
quite reasonable quantitative predictions when compared to numerically exact 
simulations of the same contact problem. 

Another physical variable of interest in contact mechanics studies is the ra- 
tio A/Afj between the real area of contact and the nominal contact area. With 
the help of Persson's theory one can derive expressions that relate A/Aq to the 
average interfacial separation u. Next, the predictions from such a relation can 
be compared to those of numerical calculations. Figure [T^] depicts the fractional 
contact ratio obtained in theoretical and numerical simulations of rough surfaces 
with variable Hurst exponents over a wide pressure region. As depicted, in the 
low-pressure regime (large u zone) all approaches converged to the same limit. 
This is further proof of the suitability of the numerical techniques discussed 
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Figure 10: Probability distribution P{u) for several squeezing pressure values. 
The red line corresponds to the theoretical predictions and the blue to numerical 
simulations (in (a) and (c) using BEM and in (b) using GFMD). 




Figure 11: Probability distribution P{u) for the squeezing pressure p = O.OOSi? 
obtained in a numerically exact calculation (blue curve), and by applying the 
theoretical result from (Eq. ([TE)) ) (red curve). Also shown is the theory predic- 
tion (Eq. PT)) ) for the distribution Pi{u) of boundary-line averaged interfacial 
separations (green curve). 
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here to study the contact mechanics of elastic solids with randomly-rough sur- 
faces under engineering conditions. As the pressure and H exponent increase, 
discrepancies arised between GFMD and the other approaches. In future works 
the cause for such discrepancies in the high pressure region must be further 
investigated. 




Figure 12: The contact area A/Aq as a function of average interfacial separation 
u for the H = 0.3, 0.5 and 0.8 surfaces. 



5. Summary and conclusions 

The contact mechanics theory developed by Persson has been extended to 
allow for the calculation of the distribution of interfacial separations P{u). The 
theory has been applied to study the contact mechanics of a flat elastic solid 
squeezed against an infinitely-hard randomly rough substrate. Three different 
coarse-grained numerical approaches have been used to simulate the same prob- 
lem: (a) the boundary element method (BEM), (b) Green's function molecular 
dynamics (GFMD) and (c) smart-block molecular dynamics. The theoretical 
predictions have been compared to those of the numerical methods. 

All the numerical methods and the analytical theory gives very similar re- 
sults for the pressure p and the fractional contact area A/Aq, as a function of 
the average interfacial separation u. In agreement with some earlier numerical 
studies when find a linear proportionality between real area of contact A and the 
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applied load at low loads, and a logarithmic relation between the average interfa- 



cial separation and the applied pressure in the same load regime (lAkarapu et aL. . 
2OIOI: ICampana et all 1201 it iLorenz fc Perssonl 120091: iLorenz et all |2010| ). We 



note that these functional forms are those predicted by Persson's theory thus 
pointing to the capabilities of the theory to properly account for long-range 
elastic deformations. 

The distribution of interfacial separations P{u) was studied in the low-to- 
medium pressure regions and our results showed that the theory gives quantita- 
tive predictions of reasonable accuracy for P{u) when compared to the results 
obtained from numerically-exact calculations. However, the rather small sys- 
tem sizes used in the numerical calculations resulted in finite size effects (noise) 
for large interfacial separations. Among the numerical schemes we note that 
the BEM model, based on a FFT-accelerated implementation of the Boussinesq 
equation, is fast and accurate over the whole range of squeezing pressures. The 
GFMD method is also computationally very fast, but its results deviated from 
the expected solution for large values of the roughness exponent and high pres- 
sures. While the cause for such discrepancies needs to be investigated, we note 
that it lies within the current numerical implementation of the GFMD code and 
not in the GFMD theory which is in principle exact. Lastly, the smart-block 
classical MD is the most computationally demanding approach. Nevertheless, 
in contrast to BEM and GFMD, it naturally includes adhesion and friction in 
an atomistic way, and its current implementation can be applied within the full 
pressure range. 
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